  //////////////////////////////////////////////////
  // MSM version 1.12
  //Multispecies Statistical Model    
  //////////////////////////////////////////////////

// In this version we add the M1 parameter
// we reduced suitabilities to include from age 1 to age 12
// two species,  this version estimates M2 
// updated 2011, "EQ" refers to corresponding equations in Jurado Molina et al 2002

DATA_SECTION
  number MIXM2
  number niter
  int iter
  !!MIXM2=.9;
  !!niter=4;
  init_int mode // If 0 single species otherwise multi
  init_int styr
  init_int nyrs
  init_int nspp
  init_int nages
  init_3darray obs_catch(1,nspp,1,nyrs,1,nages)
  init_4darray wt(1,nspp,1,nspp,1,nages,1,nages) // prey weight in the stomach of the predator
  init_3darray food_in(1,nspp,1,nyrs,1,nages)  
  init_3darray avg_growth(1,nspp,1,nyrs,1,nages)   
  init_vector other_food(1,2)  
  init_4darray suit_main(1,nspp,1,nspp,1,nages,1,nages) 
  init_matrix suit_other(1,2,1,nages)
  init_matrix M1(1,nspp,1,nages)
  init_4darray suit_transp(1,nspp,1,nspp,1,nages,1,nages)

  3darray pred_food_ration(1,nspp,1,nyrs,1,nages)  
  matrix suit_wt(1,nspp,1,nages)
  vector trend(1,nyrs)
  !! trend.fill_seqadd(-1,2/(nyrs-1));
  int k
  int i
  int j
  int nselages
  int p  // predator spp index  

  number tau
  3darray p_obs(1,nspp,1,nyrs,1,nages)
  matrix  tc_obs(1,nspp,1,nyrs)

 LOCAL_CALCS
   tau = 50.;
  
   for (k=1;k<=nspp;k++)
   {
     for (i=1;i<=nyrs;i++)
     {
       tc_obs(k,i)= sum(obs_catch(k,i));
       p_obs(k,i) = obs_catch(k,i)/(tc_obs(k,i)+.01);
       
     }
   }
 END_CALCS
   // Jim's data additions here (for surveys)
   init_ivector   n_srv_biom(1,nspp)
   init_imatrix yrs_srv_biom(1,nspp,1,n_srv_biom)
   init_matrix      srv_biom(1,nspp,1,n_srv_biom)
   init_matrix   srv_biom_se(1,nspp,1,n_srv_biom)
   matrix   srv_biom_lse(1,nspp,1,n_srv_biom)
  !! srv_biom_lse = elem_div(srv_biom_se,srv_biom);
  !! srv_biom_lse = sqrt(log(square(srv_biom_lse) + 1.));
  !! cout<<srv_biom_lse<<endl;

   init_ivector      n_srv_age(1,nspp)
   init_ivector   srv_age_type(1,nspp)
   init_imatrix    yrs_srv_age(1,nspp,1,n_srv_age)
   init_ivector   srv_age_bins(1,nspp)
   init_matrix       srv_age_n(1,nspp,1,n_srv_age)              // Sample size for multinomial
   init_3darray        srv_age(1,nspp,1,n_srv_age,1,srv_age_bins) // observed age/size compositions

   !!for (k=1;k<=nspp;k++) for (i=1;i<=n_srv_age(k);i++) srv_age(k,i) /= sum(srv_age(k,i));
   init_3darray    srv_age_err(1,nspp,1,nages,1,srv_age_bins) // observed age/size compositions
   init_number test_read;
   !! if (test_read != 12345) {cout<<"Read file error"<<endl<<test_read<<endl<<srv_age_err<<endl;exit(1);}

 INITIALIZATION_SECTION
  srv_sel_inf 4.5
  srv_sel_slp .9
  mean_rec 2.
  log_srv_q 0.;
  mean_F .2
  sel_coff 1.
  alpha 1

PARAMETER_SECTION
  init_vector mean_rec(1,nspp,1)
  init_bounded_vector mean_F(1,nspp,.001,2,1)
  init_bounded_matrix rec_dev(1,nspp,1,nyrs,-10,10,2)
  init_bounded_matrix init_dev(1,nspp,2,nages,-10,10,2)
  init_bounded_matrix F_dev(1,nspp,1,nyrs,-4,4,2)
  init_bounded_number alpha(.49,1.01,-3)
  !! nselages=8;
  init_matrix sel_coff(1,nspp,1,nselages,4)
  matrix sel(1,nspp,1,nages)
  number M2_pen;
  matrix R(1,nspp,1,nyrs)
  matrix biomass(1,nspp,1,nyrs)
  sdreport_vector depletion(1,nspp)
  3darray p_hat(1,nspp,1,nyrs,1,nages)
  matrix tc_hat(1,nspp,1,nyrs)
  init_vector log_srv_q(1,nspp,3);
  init_vector srv_sel_inf(1,nspp,-3);
  init_vector srv_sel_slp(1,nspp,-5);
  init_matrix srv_sel_coff(1,nspp,1,nselages,3)
  vector       avg_srv_sel(1,nspp);
  vector           avg_sel(1,nspp);
  vector             srv_q(1,nspp);
  matrix      srv_biom_hat(1,nspp,1,n_srv_biom)
  matrix           srv_sel(1,nspp,1,nages)
  3darray      srv_age_hat(1,nspp,1,n_srv_age,1,srv_age_bins) // observed age/size compositions
  3darray  catch_hat(1,nspp,1,nyrs,1,nages)
  3darray  residuals(1,nspp,1,nyrs,1,nages)
  3darray  N(1,nspp,1,nyrs+1,1,nages)
  3darray  avail_food(1,nspp,1,nyrs,1,nages) 
  3darray  AvgN(1,nspp,1,nyrs,1,nages)
  3darray  F(1,nspp,1,nyrs,1,nages)
  3darray  Z(1,nspp,1,nyrs,1,nages)
  3darray  S(1,nspp,1,nyrs,1,nages)
  3darray    M2(1,nspp,1,nyrs,1,nages)
  3darray M2old(1,nspp,1,nyrs,1,nages)
  objective_function_value obj_fun;
  
PROCEDURE_SECTION
  srv_q = mfexp(log_srv_q);
  M2_pen.initialize();
  M2old.initialize();
  M2.initialize();
  calc_selectivity();
  calc_numbers_at_age(1);
  if (mode)
  {
   for (iter=1;iter<=niter;iter++)
    calc_numbers_at_age(2);
  }

  for (k=1;k<=nspp;k++)
  {
    calc_predicted_values(k);
    calc_objective_function(k);
  }

FUNCTION calc_selectivity
  for (k=1;k<=nspp;k++)
  {
    avg_sel(k)               = log(mean(mfexp(sel_coff(k))));
    sel(k)(1,nselages)       = mfexp(sel_coff(k));
    sel(k)(nselages+1,nages) = sel(k,nselages);

    avg_srv_sel(k)            = log(mean(mfexp(srv_sel_coff(k))));
    if (active(srv_sel_slp)) // Asymptotic survey selectivty
    {
      for (i=1;i<=nages;i++)
        srv_sel(k,i) = 1/(1+exp(-srv_sel_slp(k)*(double(i)-srv_sel_inf(k))));  // EQ 7: Sa=1/(1+exp(u-va)) (modified)
    }
    else  //survey selectivity uses coefficients
    {
        srv_sel(k)(1,nselages)       = mfexp(srv_sel_coff(k));
        srv_sel(k)(nselages+1,nages) = srv_sel(k,nselages);
        srv_sel(k)                  /= mean(srv_sel(k) );
    }
  }

FUNCTION calc_available_food
  for (p=1; p<= nspp; p++)  // predator loop
    {
    for (j=1;j<=nages;j++)  // predator age loop
    {
      for (k=1;k<=nspp;k++)  // prey loop
      {
        for (int prey_age =1;prey_age<=nages;prey_age++)  // prey age loop
        {
          avail_food(p,i,j) += suit_main(p,k,j,prey_age)* AvgN(k,i,prey_age)* wt(p,k,j,prey_age); // EQ: 2 (dnom) M2=small((Npred*R*Suit)/(sum(Nprey*W*suit))
       } //end prey age loop
      }  // end prey spp loop
      avail_food(p,i,j) += other_food(p)* suit_other(p,j); 
    } // Pred age (j)
  }  // pred loop
 
FUNCTION calc_M2
   
  for (int prey_k=1;prey_k<=nspp;prey_k++)  // prey spp loop
  {
    for (int prey_age =1;prey_age<=nages;prey_age++)  // prey age loop
    {
      dvariable Mtmp=0.;
      for (p=1;p<=nspp;p++)   // pred species loop
      {  
        for (j=1;j<=nages;j++)  // Pred age loop
        {
          Mtmp += AvgN(p,i,j) * food_in(p,i,j)*suit_main(p,prey_k,j,prey_age)/avail_food(p,i,j); // EQ: 2 (numr) R="food_in" M2=small((Npred*R*Suit)/(sum(Nprey*W*suit))
        }  // end pred age loop
        M2(prey_k,i,prey_age) = Mtmp;
     }  // end pred spp loop
    }   // end prey age loop
  }   // end prey spp loop
   
 /* for (k=1;k<=nspp;k++)   // pred species loop
  {
    M2(k,i)    += elem_div(elem_prod(AvgN(k,i),suit_main( avail_food(k,i)); 
    M2(k,i)    = MIXM2 * M2old(k,i)  + (1.-MIXM2)*M2(k,i);
    M2old(k,i) = M2(k,i);
  }
  */ 

FUNCTION calc_mortality
  F(k,i) =  sel(k) * mean_F(k) * mfexp(F_dev(k,i)); // EQ 6 Fa=Sa*Ft
  Z(k,i) = F(k,i) + M1(k)(1,nages) + M2(k,i);
  S(k,i) = mfexp(-Z(k,i));

FUNCTION void calc_numbers_at_age(int pass_number)
  for (k=1;k<=nspp;k++)
     {
     // Sub-vector of rec_dev(k) from 1 to nyrs (rec_dev is longer now...)
     R(k) = mfexp(mean_rec(k) + rec_dev(k)(1,nyrs) ); // Recruitment
     // Top row, left column....
     for (i=1;i<=nyrs;i++)
        N(k,i,1) = R(k,i);      // Recruitment  
        // Initial age composition (first year)
        for (j=2;j<=nages;j++)
           N(k,1,j)= mfexp(mean_rec(k) + init_dev(k,j) );
    switch (pass_number)
    {
      case 1:
       AvgN(k,1) = N(k,1);
       break;
      case 2:
      default:
      AvgN(k,1) = elem_div( elem_prod( N(k,1),(1.-S(k,1))) ,Z(k,1)); 
       break;
     }
  } // End loop over species

  // Main year loop for filling numbers at age
  avail_food.initialize();
  for (i=1;i<nyrs;i++)
     {
      if (mode && current_phase()>2)
        {
         calc_available_food();
         if (pass_number>1)
         calc_M2();
        }
      for (k=1;k<=nspp;k++)
         {
          calc_mortality();
          N(k,i+1)(2,nages)= ++elem_prod(N(k,i)(1,nages-1),S(k,i)(1,nages-1)) ; // EQ 4 Na+1,t+1,q+1=Natq*exp(-(Fat+Ma))
          N(k,i+1,nages) += S(k,i,nages)*N(k,i,nages);
          switch (pass_number)
         {
       case 1:
        AvgN(k,i+1) = N(k,i+1);
        break;
       case 2:
       default:
        AvgN(k,i+1) = elem_div( elem_prod( N(k,i+1),(1.-S(k,i+1))) ,Z(k,i+1));
        break;
      }
    }
  }
  i=nyrs; // for the terminal year
  if (mode && current_phase()>2 )
  {
    calc_available_food();
    calc_M2();
  }

  for (k=1;k<=nspp;k++)
  {
    calc_mortality();
    N(k,i+1)(2,nages)= ++elem_prod(N(k,i)(1,nages-1),S(k,i)(1,nages-1)) ; // EQ 4 Na+1,t+1,q+1=Natq*exp(-(Fat+Ma))
    N(k,i+1,nages) += S(k,i,nages)*N(k,i,nages);
  }

FUNCTION void calc_predicted_values(int k)
  catch_hat(k) = elem_prod(elem_div(F(k),Z(k)) ,elem_prod(1.-mfexp(-Z(k)) , N(k))); // EQ 5 Cat=(F*N*(1-exp(-Z))/(F+M1+M2)
  for (i=1;i<=nyrs;i++)
  {
    tc_hat(k,i)= sum(catch_hat(k,i));
    p_hat(k,i) = catch_hat(k,i)/tc_hat(k,i);
    biomass(k,i) = N(k,i) * avg_growth(k,i)(1,nages);
  }
  
  for (i=1;i<=n_srv_biom(k);i++)
  {
    // convert years into indices for 1-=nyrs counting purposes 
    int yr_ind = yrs_srv_biom(k,i) - styr + 1;
    // note: assumes survey occurs in mid-year (pow(,.5))
    srv_biom_hat(k,i) = elem_prod(srv_q(k) * srv_sel(k) , elem_prod(pow(S(k,yr_ind),0.5),N(k,yr_ind)) ) * avg_growth(k,yr_ind)(1,nages);  
    // EQ 9 sort of Bhat=q*Sa*sqrt(Survival)*N*Wt
  }
  for (i=1;i<=n_srv_age(k);i++)
  {
    // convert years into indices for 1-=nyrs counting purposes 
    int yr_ind = yrs_srv_age(k,i) - styr + 1;
    // need to test for type of age data here (if length bins are different
    dvar_vector tmp_age = elem_prod(srv_sel(k) , elem_prod(pow(S(k,yr_ind),0.5),N(k,yr_ind)) );
    if ( srv_age_type(k)==1)
      srv_age_hat(k,i) = tmp_age / sum(tmp_age);
    else
    {
      srv_age_hat(k,i) = tmp_age * srv_age_err(k);
      srv_age_hat(k,i) /= sum(srv_age_hat(k,i));
    }
  }

 // OjO what's up with this?
  if (sd_phase())
  {
    depletion(k) = biomass(k,nyrs)/(1.e-20 + biomass(k,1));
    for (i=1;i<=nyrs;i++)
    {
      // prey_consumed(k,i) = wt(k,i)*elem_prod(M2(k,i),AvgN(k,i)); // Vector * vector = scalar
      //prey_consumed(k,i) = wt(k,i)(1,nages)*elem_prod(M2(k,i),AvgN(k,i)); // Vector * vector = scalar
    }
  }

FUNCTION void calc_objective_function(int k)
    //residuals(k)= elem_div(obs_catch(k)-catch_hat(k),sqrt(catch_hat(k) + .1));
    // Multinomial
    obj_fun -= tau*sum(elem_prod(p_obs(k),log(p_hat(k)+1.e-4)));

  // Errors in total catch estimation
    obj_fun += 100.*norm2(log(tc_obs(k))-log(tc_hat(k)+1e-4));

    if (!last_phase())
      obj_fun+= 50.*square(log(mean(F(k))/.2));
    else
      obj_fun+= 10.*square(log(mean(F(k))/.2));

  // Invoke a penalty when the partial F's go down with age
    for (j=1;j<=nages-1;j++)
      if (sel(k,j)>sel(k,j+1))
        obj_fun += 20.*square(log(sel(k,j)/sel(k,j+1)));

    obj_fun += 10.*square(avg_srv_sel(k)); // this part will go to zero (it's a condition)
    obj_fun += 10.*square(avg_sel(k)); // this part will go to zero (it's a condition)
    obj_fun += 10. * norm2(first_difference(first_difference(log(srv_sel(k)))));
    obj_fun += 10. * norm2(first_difference(first_difference(log(sel(k)))));

  // cout <<srv_sel<<endl<<srv_q<<endl<<srv_biom_hat<<endl;
  // for (k=1;k<=nspp;k++)
  for (i=1;i<=n_srv_biom(k);i++)
    obj_fun += square(log(srv_biom(k,i)) - log(srv_biom_hat(k,i)) ) / (2.*srv_biom_lse(k,i)*srv_biom_lse(k,i));
        
  for (i=1;i<=n_srv_age(k);i++)
    obj_fun -= srv_age_n(k,i)*sum(elem_prod(srv_age(k,i),log(srv_age_hat(k,i)+1.e-4)));
  obj_fun +=  1.* norm2(rec_dev(k));
  obj_fun +=  1.* norm2(init_dev(k));
  obj_fun +=  1.* norm2(F_dev(k));
  obj_fun += 100.*(square(mean(F_dev(k))) + square(mean(rec_dev(k))) + square(mean(init_dev(k))));
  //if (active(F_dev))
   // obj_fun += 10.* square(trend*F_dev(k))/norm2(F_dev(k)+1e-20);

  //obj_fun += 10.* square(alpha-1);
  obj_fun += 10.* M2_pen;

REPORT_SECTION
  cout <<endl<<"========End of phase: "<<current_phase()<<" ============"<<endl<<endl;

  report << "INPUTS"<<endl;
  report << "cATCH AT AGE"<<endl;
  report << obs_catch <<endl;

  report << "weight at the stomach"<<endl;
  report << wt <<endl;
 
  report << "food_in"<<endl;
  report << food_in <<endl;

  report << " avg weight"<<endl;
  report << avg_growth <<endl;
  
  report << "other food"<<endl;
  report << other_food <<endl;

  report << "suit plk as predator, plk as prey, age of prey 1?"<<endl;
  report << suit_main(1) <<endl;
  report << "suit cod as predator, plk as prey, age of prey 1?" << endl;
  report << suit_main(2) <<endl;

  report << "suitabilities other food"<<endl;
  report << suit_other <<endl;

  report << "residual mortality"<<endl;
  report << M1 <<endl;  

  report << "suit transp plk as predator, plk as prey, age of prey 1?"<<endl;
  report << suit_transp(1) <<endl;
  report << "suit transp cod as predator, plk as prey, age of prey 1?" << endl;
  report << suit_transp(2) <<endl;
  //report << "food intake" << endl;
  //report << food_in(1) <<endl;
  //report << "food intake" << endl;
  //report << food_in(2) <<endl;

  report << "OUTPUTS"<<endl;

  report << "Numbers at age" <<endl;
  report << N<<endl<<endl;
  report << "Fishing mortality at age" <<endl;
  report << F<<endl<<endl;
  report << "Survival at age" <<endl;
  report << S<<endl<<endl;

  //report << "Residuals"<<endl;
  //report << residuals  <<endl<<endl;
  report << "Biomass"<<endl<<biomass<<endl;

  report << "Predicted total catch"<<endl;
  report <<  tc_hat <<endl;
  report << "observed total catch"<<endl;
  report <<  tc_obs <<endl;

  report << "Predicted catch at age"<<endl;
  report <<  p_hat <<endl;
  report << "observed catch at age"<<endl;
  report <<  p_obs <<endl;
  report << "Values for M2        "<<endl;
  report << "plk          "<<endl;
  report <<  M2(1) <<endl<<endl;


  report << "pcod          "<<endl;
  report <<  M2(2) <<endl<<endl;


  report << "available food plk"<<endl;
  report << avail_food(1)<< endl;
  report << "available food cod"<<endl;
  report << avail_food(2)<< endl;
  report << "average N plk"<<endl;
  report << AvgN(1)<< endl;
  report << "average N cod"<<endl;
  report << AvgN(2)<< endl;
  report << "======Survey_selectivity-at-age==========="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    report << srv_sel(k)<<" "<<endl;
  }
  report << "======Survey_Biomass_Fit=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=n_srv_biom(k);i++)
      report << yrs_srv_biom(k,i)<<" "<<srv_biom(k,i)<< " "<<srv_biom_hat(k,i)<<endl;
  }
  report << "======Survey_age_composition_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=n_srv_biom(k);i++)
      report << yrs_srv_age(k,i)<<" "<<srv_age(k,i)<< " "<<srv_age_hat(k,i)<<endl;
  }
  report << "======Fishery_age_composition_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    for (i=1;i<=nyrs;i++)
      report << styr+i-1 <<" "<<p_obs(k,i)<< " "<<p_hat(k,i)<<endl;
  }
  report << "======Total_catch_fits=================="<<endl;
  for (k=1;k<=nspp;k++)
  {
    report << "==Species: "<<k<<"----"<<endl;
    report << tc_obs(k)<< endl<<tc_hat(k)<<endl;
  }

RUNTIME_SECTION
  convergence_criteria .1, 0.001,0.00000001
  maximum_function_evaluations 200, 300, 3000

TOP_OF_MAIN_SECTION
  gradient_structure::set_MAX_NVAR_OFFSET(1000);
  gradient_structure::set_GRADSTACK_BUFFER_SIZE(3000000);
  gradient_structure::set_CMPDIF_BUFFER_SIZE(4000000);






